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FINAL  REPORT 


Project  Identification:  AFOSR  Grant  F49620-96- 1-0133 

Project  Title:  Statistical  Modelling  and  Simulation  for  Microstructure  in  Materials  Science 
Principal  Investigator:  Chuanshu  Ji,  UNC-Chapel  Hill,  cji@stat.unc.edu 

The  main  accomplishment  in  this  project  consists  of  two  parts:  (I)  Reconstruction  cycle  in  mi¬ 
crostructure  modelling;  and  (II)  Computation  of  effective  properties  via  Markov  chain  Monte  Carlo 
(MCMC)  algorithms. 

The  details  are  contained  in  Derr  (1998)  —  a  draft  of  Ph.D.  dissertation  entitled  “Statistical 
modelling  of  microstructure  with  applications  to  effective  property  computation  in  materials  science” 
by  Bob  Derr  under  the  supervision  of  Chuanshu  Ji.  Several  papers  based  on  this  thesis  are  being 
written  and  to  be  sumitted  for  publication. 

I.  Reconstruction  cycle  in  microstructure  modelling 

One  of  the  most  important  issues  in  materials  science  is  the  connection  between  materials  prop¬ 
erties,  e.g.  conductivity,  elastic  moduli,  strength,  etc.  and  microstructures.  Along  this  line,  many 
computer  models  were  proposed  to  generate  synthetic  microstructures  on  which  some  numerical 
schemes,  e.g.  finite  element,  were  used  to  calculate  the  materials  properties  of  interest,  assuming  the 
local  properties  satisfy  certain  partial  differential  equations  for  stress/strain  or  diffusions.  On  the 
other  hand,  experimentation  was  conducted  in  laboratories  which  measures  those  properties  from 
real  materials.  A  significant  gap  exists  between  these  two  aspects  of  the  study  due  to  the  lack  of 
methodology  for  fitting  the  computer  models  (i.e.  estimating  the  parameters  in  those  models)  based 
on  the  real  microstructure  data. 

Several  steps  were  taken  in  this  project  towards  bridging  such  a  gap.  A  class  of  statistical  models 
for  microstructure  —  hard-core/soft-shell  elliptical  processes,  simply  called  elliptical  processes,  was 
proposed,  and  the  related  modelling  task  was  carried  out  via  the  following  reconstruction  cycle: 
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Step  1.  Descriptive  data  summary:  microstructure  image  processing; 


Step  2.  Model  fitting:  parameter  estimation  in  the  proposed  probability  density; 

Step  S.  Simulation  from  the  fitted  model:  Markov  chain  Monte  Carlo  (MCMC)  algorithms; 

Step  4.  Model  checking:  good  of  fit  tests,  comparisons  of  features. 

Although  this  follows  from  the  basic  principle  in  statistical  modelling,  its  implementation  on 
models  as  complex  as  the  elliptical  processes  is  new  and  challenging.  The  main  consideration  in 
choosing  this  model  class  is  the  combination  of  practicality  and  feasibility.  The  hard-core  assumption 
is  essential  in  modelling  microstructure  since  particles  (grains)  in  3D  are  nonoverlapping.  The 
main  technical  difficulty  with  the  hard-core  assumption  lies  in  that  it  induces  complicated  spatial 
dependence  structures.  The  class  of  elliptical  processes  has  several  advantages.  On  one  hand,  it 
is  rich  enough  to  generate  a  great  variety  of  synthetic  microstructural  patterns  (see  Fig.  lc,  Id, 
2a  -  2f),  because  it  allows  variability  in  location,  orientation,  aspect  ratio  and  size  of  ellipses. 
On  the  other  hand,  it  is  computationally  manageable  (although  intensive),  since  it  is  still  a  finite¬ 
dimensional  model  (the  dimensionality  can  be  quite  high).  More  realistic  microstructures  and  flexible 
deformable  shapes  can  be  approximated  when  computational  facilities  permit.  Although  this  report 
only  summarizes  results  for  the  2D  elliptical  processes,  the  extension  to  3D  ellipsoidal  processes  can 
be  handled  in  a  similar  manner. 

Step  1.  Take  the  real  microstructures  in  Fig.  la  and  lb,  approximate  them  by  nonoverlapping 
ellipses  via  range-oriented  fitting  or  least-squares  fitting  [cf.  Derr  (1998)]. 

Step  2.  A  configuration  x  —  (xi , . . .  ,  xn}  consists  n  ellipses,  each  denoted  by  xt  =  (£«,  &,-o,i-bt), 
with  center  &,  orientation  0,,  major  and  minor  half-axes  Oj  and  respectively.  The  probability 
distribution  of  an  elliptical  process  is  characterized  by  its  density  f(x)  with  respect  to  a  reference 
measure  p  x  v.  Under  p,  {£»}  form  a  homogeneous  Poisson  process  in  a  bounded  region  S  C  M2 
with  the  mean  measure  being  the  area  m(S)  of  5,  and  under  v,  { [9i ,  a, ,  6t) }  are  iid  marks  with  0* 
uniformly  distributed  on  [0, 7r],  and  (oi,  bt)  uniformly  distributed  over 'the  triangle  {l  <  bi  <  a*  <  L}, 
where  l  <  L  are  two  known  positive  constants,  to  be  specified  in  the  experiment. 

The  density  f(x)  has  a  form 

/(x;  A,  a,  b)  =  A  n/i(x;  a,  b)[B{\,  a,  6)]-1, 
with 

(1) 

h{x;  a,  b)  =  I(b<bni<ann<a)Mx)Y[9(xi)i 

i=l 


where 

(i)  ann  =  max  dj,  and  6ni  =  min  bi  [in  general,  anj  denotes  the  ith  order  statistics  in  (01, .  ■ .  ,a„)], 

l<i<n  l<t<n 

and  I  a  is  the  indicator  function  of  an  event  A ; 

(ii)  A,  a  and  b  are  three  unknown  parameters,  with  the  assumed  ranges  A  >  0  and  l  <  b  <  o  <  L, 
and  the  sufficient  statistics  {n,Onn,6n  1}; 
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(iii)  g{xi)  >  0  Vx*  and  Vi,  with  a  known  function  g,  depends  on  Xi  only  through  its  area,  e.g. 
g(xi)  =  1,  or  g(xi)  =  Oj&i,  or  g(xi)  =  (a^i)-1,  to  give  a  few  particular  choices; 

(iv)  The  factor  A(x)  >  0  Vx  includes  the  range  restriction  0  <  8,  <  ir,  a,  >  b,  Vi,  and  the 
nonovelapping  restriction,  etc.,  but  does  not  involve  the  unknown  parameters; 

(v)  The  normalization  factor  is  expressed  as 

B(X,a,b) 

(2)  =1 l{e-m^S)[m(S)]k/k\}  f  [  f  Afe/l(x;a,6)f[(deid0idaid6<), 

k  JSk  J[0,*)k  J[b,a]^  tti 

with  the  term  e~m(~s )  when  k  =  0.  Note  that  B(X,  a,  b)  is  practically  a  finite  sum.  Since  every 
ellipse  has  a  minor  half-axis  no  less  than  l  >  0,  and  the  ellipses  are  nonoverlapping  with  their 
centers  distributed  over  the  bounded  region  S,  any  configuration  consists  only  a  finite  number  of 
such  ellipses. 

Denote  the  maximum  likelihood  estimates  (MLE)  for  A,  a  and  b  by  A,  a  and  b  respectively.  Then 

(3)  a  =  ann  and  S  =  6n  i> 


The  MLE  A  can  be  found  by  solving  the  equation 


(4) 


A  S(A,  onn,  6ni)  nB(X,ann,bni)  —  0, 


using  Monte  Carlo  approximate  likelihood  [cf.  Derr  (1998)].  Alternatively,  the  ratio  appears 
to  be  a  simple  reasonable  estimate  of  A. 


Step  3.  To  simulate  spatial  patterns  from  a  target  distribution  defined  on  a  large  configuration  space 
Cl,  the  basic  idea  of  MCMC  is  to  construct  a  Markov  process  X  =  {Xt}  such  that  it  has  the  target 
distribution  as  its  (unique)  equilibrium  distribution.  The  key  is  to  incorporate  the  knowledge  of  the 
target  distribution  into  the  transition  probability  mechanism  of  X.  A  special  case  of  MCMC  — 
spatial  birth-and-death  processes  (SBDP)  —  is  suitable  for  simulation  of  the  elliptic  processes.  A 
sufficient  (reversibility)  condition  to  guarantee  the  convergence  to  the  equilibrium  density  /,  with  a 
birth  rate  &(•)  and  a  death  rate  D(-),  is  the  following  detailed  balance  equation 


(5)  /(xn)  b(xn,xn+ 1)  =  f(xn+1)  D(xn,xn+i)  whenever  /(sn+1)  >  0, 


where  the  notation  xn  =  {xi, ...  ,xn}  indicates  the  number  of  objects  in  x,  and  the  dependence  of 
/  on  parameters  [as  shown  in  (1)]  is  suppressed. 

Even  with  the  knowledge  of  a  fully  specified  /  in  (5),  there  are  still  many  different  choices  for  b(-) 
and  D(-).  One  is  to  specify 


D(xn,xn+i)  =  l/5(xn+i), 
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(6) 


thus  the  birth  rate  accordingly  by  (5) 


(7) 


b(xn,xn+1)  =  D(xn,xn+1)  f(xn+1)/f(xn) 


{A,  if  xn+i  does  not  overlap  Xi  V*  =  1, . . .  ,  n, 
0,  otherwise. 


This  can  be  easily  implemented. 

Step  4-  Checking  whether  data  fit  a  proposed  model  can  be  done  either  in  the  parameter  space  or  in 
the  sample  space.  For  the  problem  considered  in  this  work,  the  MLEs  were  derived  which  enjoy  many 
good  asymptotic  properties;  but  comparisons  in  the  sample  space  might  be  more  meaningful  since 
one  really  wants  to  “see”  if  the  synthetic  microstructures  resemble  the  real  ones.  The  sample  space 
of  microstructural  images  has  high  dimensionality.  Some  average-based  norms  such  as  L2-norm,  etc. 
often  fail  to  highlight  the  difference  in  certain  particular  features  which  may  be  more  relevant  to 
the  underlying  materials.  Hausdorff  metric  for  shape  comparison  could  be  used  but  instead,  some 
simple  feature  extraction  was  used  in  this  project.  For  each  of  the  features  (location,  orientation, 
aspect  ratio  and  size),  two  histograms  were  plotted,  one  for  the  original  microstructure  image,  one 
for  the  synthetic  one  generated  from  the  fitted  model  [cf.  Derr  (1998)].  Some  goodness  of  fit  tests 
were  conducted  to  compare  the  empirical  distributions. 

Notice  that  there  are  many  other  mothods  in  materials  science  literature  to  generate  synthetic 
microstructures,  but  they  suffer  from  a  common  drawback  —  lack  of  reconstruction  capability, 
because  they  are  not  tuned  to  real  microstructures.  The  approach  taken  in  this  project  is  expected 
to  make  a  notable  impact  to  this  area. 


II.  Computation  of  effective  conductivity  via  Markov  chain  Monte  Carlo 

The  great  variability  of  local  properties  in  heterogeneous  materials  makes  the  computation  of 
global  properties  an  important  and  challenging  project.  The  result  obtained  in  this  project  concerns 
computation  of  effective  conductivity  using  MCMC.  The  approach  looks  promising  and  should  work 
for  computation  of  effective  elastic  moduli  as  well. 

The  approach  was  originally  introduced  by  in  Brown  (1955)  and  then  extended  in  a  number 
of  papers  by  Torquato  et  al.  This  approach  demonstrates  explicitly  the  connection  between  the 
effective  conductivity  and  the  statistical  distribution  of  microstructural  features,  in  particular  the 
fc-point  probability  functions  (k  =  1,2,...),  hence  it  is  capable  of  handling  multiple  representations 
in  different  length  scales.  It  becomes  quite  clear  that  the  promise  of  this  approach  relies  on  some 
ingenious  computational  techniques. 

For  the  simple  illustration  purpose,  consider  a  two-phase  composite  material  represented  by  the 
2D  domain  S  and  follow  the  notation  in  Part  I.  The  phase  i  itself  is  homogeneous  with  constant 
electrical  conductivity  a*  >  0,  i  =  1, 2.  Assume  the  phase  1  region,  J?j  C  S,  is  precisely  covered 
by  the  configuration  x  with  n  nonoverlapping  grains  (e.g.  ellipses),  and  R2  =  S\R\  is  the  phase  2 
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region.  Define  the  local  conductivity 


(8)  cr{z)  =  <7j{zeRl j  +  a2I{zeR3},  z  e  S. 

Hence  cr(z),  z  &  S  is  a  random  field  resulting  from  the  random  distribution  of  x. 

Materials  scientists  axe  interested  in  the  effective  behavior  of  heterogeneous  media.  In  this  setting, 
if  there  is  a  homogeneous  medium,  the  effective  medium  S* ,  whose  conductivity  er*  is  close  to  the 
conductivity  of  S  when  measured  on  a  large  space  scale,  then  a *  is  called  the  effective  conductivity. 

There  are  two  issues.  The  first  one  is  to  justify  the  exitence  of  a *  —  the  averaging,  called 
homogenization,  takes  place  so  that  the  complex  small  scale  random  structure  is  replaced  by  an 
asymptotically  equivalent  homogeneous  deterministic  structure.  Papanicolaou  (1995)  is  an  excellent 
survey  for  the  related  homogenization  theory. 

The  second  issue,  computation  of  a*,  is  practically  more  relevant  in  materials  science.  Fol¬ 
lowing  the  perturbation  expansion  approach  originally  developed  by  Brown  (1955),  Torquato  and 
his  collaborators  made  significant  contributions  to  this  problem  [cf.  Torquato  (1985),  Sen  and 
Torquato  (1989),  Torquato  and  Sen  (1990),  and  Torquato  (1991)].  The  perturbation  expansion  for 
d-dimensional  media  (in  particular,  d  =  2  or  3)  given  in  these  papers  is  expressed  as 

OO 

(11)  ~~  1[cr*  +  (d  —  l)(jjU]  =  qiffij  U  —  ^  i  ^  i,j  =  1,2, 

k= 2 

where  the  tensor  coefficients  axe  represented  as  integrals: 
for  k  =  2, 

(12)  4°  =  Sjjnj  l 

where  the  subscript  6  on  the  integral  indicates  that  it  is  carried  out  with  the  exclusion  of  an  infini¬ 
tesimal  neighborhood  of  z\  (i.e.  a  tiny  d-dimensional  sphere  centered  at  z\)\ 
and  for  k  >  3, 

(13)  A*11  =  (-l)V-*  t)]  j  ■■  f  (Ti2"-Tk-i,k-Di)i*r--dzk. 

The  notation  in  (11),  (12)  and  (13): 

•  qt  is  the  volume  (area)  fraction  of  phase  i; 

•  ft  -  = _ £i  ZZi _ • 

•  Each  zm  €  5,  as  a  d-dimensional  column  vector,  has  norm  |zm|,  m  =  1, . . .  ,  k\ 

•  Tm_x,m  =  | zm  —  zm—  1 1  2 [  d  (zm  —  zm—i)(zm  —  Zm—iY  —  |zm  —  zm— 1|  U  ]  is  the  dipole- 
dipole  interaction  tensor,  where  zl  represents  the  transpose  of  z  (thus  a  row  vector),  and  U  the  unit 
dyadic  [see  Nemat-Nasser  and  Hori  (1993)  Section  15  for  the  basic  tensor  notation  and  terms],  (the 

integrand  T\2  •  •  •  Tk-\,k  •  Dlk  is  a  function  of  z\, . . .  ,  zk); 
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For  the  determinant 


(14) 
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its  entry  P}ni  m2  is  the  probability  that  all  (m2  —  mi  +  1)  points  zmi ,  zm, +1 , . . .  ,  zm2  are  contained 
in  the  phase  i  region  Ri,  where  mi,  m2  =  1, . . .  ,  fc  and  mi  <  m2. 

Numerical  computation  of  is  crucial  and  intensive,  especially  for  large  fc.  Most  of  the  previous 
works  were  limited  to  the  derivation  of  lower-order  bounds  (for  k  <  5).  These  bounds  provided  useful 
approximations  of  the  conductivity  for  certain  special  cases  of  microstructure,  but  the  result  in  the 
related  error  estimate  was  lacking  and  the  gap  between  the  available  bounds  and  the  true  value  of 
conductivity  remained  large  in  many  situations.  Therefore,  the  need  for  solving  the  long-standing 
open  problem  of  deriving  higher-order  approximations  (fc  >  5)  of  for  a  variety  of  microstructures 
became  more  urgent.  In  this  regard,  this  work  has  made  good  progress  in: 

(i)  estimating  the  fc-point  probability  function  P[  k  for  a  wide  range  of  zi, . . .  ,  z/t,  based  on  MCMC 
simulation; 

(ii)  developing  some  recursive  formulas  to  reduce  the  fc-point  probability  function  to  a  sequence  of 
two-point  probability  functions,  thus  the  computational  speed  is  significantly  improved  [cf.  Derr 
(1998)]; 

(iii)  calculating  the  high-dimensional  multiple  integral  involved  in  by  Monte  Carlo  methods. 


Plasma-sprayed  coating 


Plasma-sprayed  coating  (see  Fig.  3a  and  3b)  is  an  important  application  in  which  the  relationship 
between  microstructure  and  properties  can  be  studied.  The  approach  developed  in  this  project 
appears  promising  for  the  quatitative  study  of  coating  microstructures  and  computation  of  thermal 
conductivity.  Only  partial  results  are  available  right  now.  This  is  still  a  part  of  ongoing  research. 

In  aircraft  engines  plasma-sprayed  metallic  coatings  protect  turbine  blades  from  highly  corrosive 
environments,  and  plasma-sprayed  ceramics  insulate  other  engine  parts  from  high  temperatures  [cf. 
Herman  (1988)].  The  powdered  coating  material,  carried  in  a  stream  of  gas  such  as  argon,  is  injected 
into  the  jet  of  plasma  from  a  plasma  gun.  The  flame  accelerates  the  particles,  and  they  are  melted  by 
its  high  temperature.  The  molten  droplets  are  propelled  onto  the  target  surface,  where  they  solidify 
and  accumulate  to  form  a  thick,  tenaciously  bonded  protective  coating.  Particles  continue  to  rain 
down  at  a  rate  that  depends  on  the  area  to  be  covered  and  how  fast  the  gun  moves  over  the  surface. 
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Transmission  electron  microscopy  helps  examine  the  internal  structure  of  particles  and  coatings; 
scanning  electron  microscopy  reveals  overall  shapes  and  textures.  The  internal  structure  —  a  mosaic 
of  crystalline  grains  —  contains  many  flaws,  suggesting  that  each  particle  solidifies  extremely  quickly, 
in  perhaps  a  millionth  of  a  second.  Ceramic  coatings  in  particular  reveal  a  multitude  of  flaws.  They 
are  riddled  with  cracks  formed  as  the  ceramic  cooled  and  are  honeycombed  with  voids  filled  with  air 
trapped  in  the  deposit.  Such  flaws  can  doom  a  coating  exposed  to  mechanical  stress.  If  they  extend 
all  the  way  through  the  coating,  they  also  make  it  useless  for  protection  against  corrosion. 

Strangely,  it  is  porosity  that  suits  plasma-sprayed  ceramic  coatings  to  one  of  their  most  important 
applications:  as  thermal-barrier  coatings,  which  are  insulating  coatings  for  metal  parts  exposed  to 
very  high  temperatures  in  gas  turbines  and  other  kinds  of  engines.  For  one  thing,  porosity  increases  a 
ceramic’s  insulating  ability.  Moreover,  porosity  and  microcracks  could  enhance  strain  tolerance.  This 
is  because  a  ceramic  is  brittle  in  the  first  place,  pores  do  not  weaken  the  material  but  instead  toughen 
it  by  interrupting  the  propagation  of  the  cracks  that  inevitably  form  as  the  material  is  strained.  By 
retaining  lubricants,  pores  also  benefit  some  plasma-sprayed  coatings  designed  to  protect  against 
wear.  There  is  trade-off,  however,  that  such  voids  could  be  fatal  to  many  wear-resistamt  metallic 
coatings  and  coatings  meant  to  prevent  their  substrate  from  corroding  or  oxidizing.  Some  multi¬ 
layered  coatings  are  used  now  in  which  vacuum  plasma  spraying  yields  a  dense,  void-free  metallic 
coating,  followed  by  ceramic  thermal-barrier  coatings. 

To  control  the  degree  of  porosity  in  different  coatings  when  processing  the  materials,  an  important 
step  is  the  quantification  of  porosity.  Fig.  3a  and  3b  show  the  complex  porosity  distribution.  The 
empirical  estimate  of  porosity  was  obtained  via  image  processing.  Extensive  lab  experiments  were 
conducted  at  NIST.  However,  modelling  the  coating  microstructure  appears  very  difficult.  A  starting 
point  is  to  use  a  similar  models  to  Fig.  2c  to  simulate  the  fibre  process  in  which  nearly  horizontal 
fibres  imitate  gaps  between  coating  particles,  while  nearly  vertical  fibres  may  represent  microcracks. 
Then  the  horizontal  and  vertical  components  of  the  effective  thermal  conductivity  tensor  can  be 
calculated  using  the  method  in  Part  II.  The  numerical  calculation  can  be  compared  with  the  results 
of  lab  tests  in  the  study  of  such  coatings  as  a  thermal  barrier.  Similar  studies  in  calculation  of 
mechanical  properties  of  coatings  will  be  pursued  also. 
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Figure  1:  (a)  (Processed)  SEM  image  of  in-situ  toughened  Si3N4  ceramic  crystals 
grown  in  a  sintering  process,  (b)  Mortar  (concrete)  —  sands  mixed  in  a  cement 
paste:  the  large  irregular  grains  are  sands;  the  smaller  circular  particles  are  bubbles 
or  pores;  the  laminate  is  a  concrete  mixture;  and  there  is  a  narrow  interfacial  zone 
between  the  rocks  which  keeps  two  rocks  from  actually  touching,  (c)  Synthetic 
microstructure  simulated  from  an  elliptical  process  fitted  by  the  Si$N4  image  data 
(a),  (d)  Synthetic  microstructure  simulated  from  a  mordified  elliptical  process 
fitted  by  the  mortar  image  data  (b).  The  modification  involves  coding  all  ellipses 
into  two  different  types,  in  which  the  sands  are  outlined  while  the  pores  are  filled 
in  black. 
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Figure  2:  Tlie  fibres  (line  segments)  in  (a),  (b)  and  (c)  are  actually  narrow  ellipses 
(aspect  ratio  0.005-  with  lengths  uniformly  distributed  in  [0.01,0.1]:  (a)  randomly 
oriented  fibres;  (b)  fibre  orientation  constrained  to  be  within  tt/8  of  the  vertical; 
(c)  fibre  orientation  constrained  within  tt/16  of  the  vertical  or  horizontal.  The 
ellipses  in  (d),  (e)  and  (f)  are  samples  of  SBDP  using  different  forms  of  g(-)  in  the 
death  rale  (6):  (d)  constant  #(■);  (e)  g{ •)  inversely  proportional  to  the  size  of  the 
ellipses:  (f)  c/(-)  proportional  to  the  size  of  the  ellipses. 
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Figure  3:  (a)  is  taken  from  Geary  (1980)  which  shows  the  microstructure  of  an 
aluminum- oxide  coating  (the  upper  thick  dark  layer)  over  a  nickel- aluminum  bond 
coating  (the  lower  thin  light  layer),  (b)  is  a  microstructure  (with  width  50  fim)  of 
Zirconia  ZrO-z  coating. 
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